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NOTATIONS 


concentration of dispersing mass in liquid phase i„e. mass 
of pollutant per unit volume of solution (mg/l) 
concentration of source fluid (mg/l) 

9 

hydrodynamic dispersion tensor (m~/ m i. n ) 

2 

ongitudinal dispersion coefficient (m /min) 

2 

coefficient of molecular diffusion (m /min) 

2 

coefficient of turbulent diffusion (m /min) 

2 

transverse dispersion coefficient (m /min) 
porosity of porous media i.e. aquifer 

volumetric fluid injection rate at the source per unit 
volume of medium (min ^) 
retardation factor 

the mass of pollutant adsorbed per unit volume of media 
(mg/l) 

the time (min) 

the component of Darcy' s velocity vector in i direction 
(m/min) 

longitudinal dispersivity (m) 

transverse dispersivity (m) 

flux at seepage boundary (mg/l per meter) 

radioactive decay coefficient equal to reciprocal of mean 
life time of radioactive pollutant (min ) 
values of nodal concentration 


width of domain. 



ABSTRACT 


In the present study an attempt has been made to study the 
pollutant transport in groundwater flowing through aquifers. A 
suitable mathematical model has been developed to predict the 
spread of pollutant in groundwater. This model incorporates many 
different methods of pollutant transport like mechanical dispersion 
molecular diffusion, hydrodynamic dispersion etc. and other pheno- 
mena such as adsorption, radioactive decay etc. 

The resulting differential equations describing pollutant 
transport have been solved for various boundary conditions and 
other parameters in both - one and two dimensional, cases using 
finite element method. It is seen that numerical solution and 
the analytical results match closely for some special case. It is 
also seen that for radioactive elements with higher half life 
period the effect of radioactive decay on pollutant concentration 
is not significant. Also the volumetric fluid injection rate and 
its pollutant concentration affects the concentration profile even 
at a distance far away from the source. For two dimensional cases 
the dispersion equation has been solved for two different sets 
of boundary conditions, and the concentration contours for 
one set of boundary conditions are presented. 



CHAPTER 1 


INTRODUCTION 

1. 1. General : 

Water resources engineer is often concerned with proper 
management and effective use of groundwater. It is important for 
such a person to be able to predict the changes in a system due 
to any proposed policy and thus obtain the new state of the 
system, given the initial system. Such a prediction before hand 
helps in checking the feasibility of the system and also having 
knowledge before hand about any ill effects that it may cause. 

Also such prediction helps in choosing from several options, the 
best policy according to some criteria. For example in a ground- 
water system it should be possible to predict water levels, sali- 
nity, spring discharges, land subsidences etc., vhich are the 
system state variables which may change due to any proposed policy 
like pumping, artificial recharge etc. 

In many regions groundwater does not pose major biological 
problem or physical quality problem, but in certain cases the 
presence of faulty sewage pipes, oil spillage etc. may cause 
severe pollution of groundwater. Pollutants, thus introduced 
into gro un dwater, may travel large distances in an aquifer and 
thus pollute groundwater over a large area. Also groundwater is 
liable to pick up minerals viiile passing throucfi rocks, containing 
such minerals. Ihis may lead to deterioration of groundwater 
quality. Another problem is that of salinity. Most groundwaters 
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becomes salinised by pollutants broughtdown from the surface. In 
natural condition an equilibrium is reached between water leaving 
the formation and pollutants carried with it. But due to activi- 
ties like groundwater recharge/ water of inferior quality is 
introduced/ then the pollutant concentration in the groundwater 
increases and this leads to increase in salinity. 

1. 2. Objective of Study : 

Groundwater quality problem in water resources engineering 
involves the mass transport in porous media where 'mass 1 is some 
pollutant moving with water in the voids of soil both in saturated 
and unsaturated zones. The mechanism which affect the transport 
of pollutant are convection/ mechanical dispersion, molecular 
diffusion, solute pollutant interaction, various chemical reac- 
tions and radioactive decay phenomenon. 

In the present work the laws governing the movement and 
accumulation of pollutants in groundwater flow have been studied. 
The movement of pollutant has been described using partial diff- 
erential equations subjected to ini rial and boundary conditions. 
These equations have then been solved by finite element method 
which is a very powerful technique. The solution of such equa- 
tions enable the engineer to predict the distribution of pollu- 
tants in tiie aquifer. 

1.3. Significance of Study : 

A major problem, which is of interest in development and 
management of water resources system, is that of water quality. 
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In fact, with increased demand for water in most parts of world 
and intensification of water utilisation, the quality problem 
becomes the limiting factor in development of water resources in 
many parts of vorld. Although in such regions, the quality of 
both surface and groundwater resources deteriorates as a result 
of pollution, special attention should be devoted to the pollu- 
tion of groundwater in an aquifer due to their very slow velocity. 
Hence, although it seems that groundwater is more protected than 
surface water against pollution, it is still subjected to pollu- 
tion and when latter occurs, the restoration to ori ginal , non- 
pollutated state, is more difficult (Bear, 1972). 

1.4. Scope of Study : 

In the present work the aquifer considered is saturated, 
porous medium. The effects of convection, hydrodynamic dispersion, 
adsorption and radioactive decay on pollutants are considered in 
one dimensional and two dimensional flow fields. In one dimen- 
sional case linear and quadratic elements have been taken. Solu- 
tions of two dimensional dispersion equation with different boun- 
dary conditions have been presented. Programs have been developed 
for F.E.M. solution of one (for linear and quadratic cases) and 
two dimensional cases. 

1.5. Organisation of Report : 

The study is reported in five chapters. The sequence of 
the topics discussed along with a brief description of each 
chapter is given below. 
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In Chapter 2 the wark done previously in this field is 
reported. Ihis discussion mainly consists of a chronological 
report of the various attempts that were made at mathematical 
modeling of pollutant transport in groundwater using finite ele- 
ment method of solution. 

In Chapter 3 a brief qualitative review of the various 
processes of pollutant transport in groundwater has been given, 
followed by the governing partial differential equation which 
describes the transport both in one dimensional and two dimen- 
sional flow fields. Then the finite element formulation of the 
above equations have been given. 

In Chapter 4 the results obtained by solving the partial 
differential equations subjected bo various boundary conditions 
by the finite element method, have been described in detail . In 
one dimensional case results obtained by the finite element 
method have been compared with the results obtained analytically 
in one special case. Two dimensional problem has also been 
solved. The results have been presented in tabular and graphical 
forms. 

In the last chapter a brief summary of the work has been 
given with suggestions for future work. 
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CHAPTER 2 

LITERATURE REVIEW 

In the field of groundwater management, the prediction 
of groundwater quality by using mathematical model warrants 
special attention. Several attempts have been made in the past 
for developing suitable mathematical models to predict pollutant 
transport in groundwater. A brief review of the work done in 
this field is presented in the following paragraphs. 

Pinder (1973) had simulated the movement of groundwater 
by Galerkin method of approximation in connection with finite 
element method of analysis. In solving the groundwater flow and 
mass transport equations through this approach, he had used func- 
tional representation of dispersion tensor, transmissivity tensor 
and fluid velocity as well as an accurate representation of boun- 
daries of irregular geometry. 

Nishi et. al . (1976) considered the hydrodynamic disper- 
sion in seepage from a triangular channel. In this they employed 
isoparametric elements to solve the two dimensional convection- 
dispersion problem. He made a comparison between finite element 
and finite difference solutions to this problem. 

Narsimhan et.al. (1978) developed mixed implicit- explicit 
Galerkin finite element method vhich was ideally suited for a 
wide class of problems arising in subsurface hydrology. These 
problem include confined saturated flow, unconfined flow under 
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free conditions subjected to Dupuit* s assumption and. flow in aqui- 
fers which are partly confined and partly unconfined. 

Lindorff (1979) studied the severity and extent of ground- 
water contamination. According to him the severity and extent 
of groundwater contamination is determined by ( 1) hydrogeologic 
setting, (2) nature of contaminant and (3) effectiveness of regu- 
latory action. 

Grove and Wood (1979) tried to predict and verify the sub- 
surface water quality changes during artificial recharge. 

Noorishad and Mehran (198 2) obtained the two dimensional 
transient dispersive- convective transport of nonconservative 
solute species in porous media. They developed a tvo nodal point, 
one dimensional, transport element for fractures -vAiich provides a 
number of advantages relative to conventional fracture represen- 
tation of two dimensional continuum elements. 

Cooley (1983) presented some new procedures for forming 
numerical solutions of saturated flow problems. The techniques 
for addressing the problems, controlling the stability of non- 
linear equation solver, devising a reliable, yet efficient method 
for determining the position of seepage surface, have been applied 
in subdomain of finite element discretization of the governing 
flow equations. He presented the solutions of problems like 
(a) steady state flow through square embankment, (b) steady state 
flow to a well, (c) drainage from a square block, (d) drainage 
involving multiple seepage faces, (e) one dimensional vertical 
infiltration, (f) transient seepage from, stream and (g) steady 
state groundwater flow around lakes. 
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Huykom et.al. (1984) developed a Galerkin finite element 
formulation for the numerical solution of water flow in saturated 
soil systems. He included a solution strategy based on Picard 
and Newton algorithm. Both algorithms are formulated for both 
rectangular and triangular elements. In this a comparative study 
of Picard and Newton-Raphson algorithm is provided. Four examples 
are presented to demonstrate the effectiveness of the finite 
element approach. These examples show that the nonlinear solution 
schemes are capable of accommodating cases involving larger vari- 
ations in saturated hydraulic conductivity as well as highly 
nonlinear soil moisture characteristics. 

Yang (1988) conducted a two dimensional saturated- 

unsaturated water and solute transport experiment in a soil tank 
infiltration and drainage condition. In this, a model accounting 
for mobile and immobile water phases is used to describe the 
solute transport. A Galerkin finite element model is designed 
to solve the transient movement of water in s a turated-un saturated 
soil. The seepage velocities at nodal points are obtained by 
using the same finite element as for the pressure field. 
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CHAPTER 3 

MATHEMATICAL PROBLEM FORMULATION 


3.1. General : 

Pollutant is a term used to denote dissolved matter 
carried with the groundwater and accumulating in the aquifer. 

For groundwater pollution study, water resources engineer must 
be familiar with potential sources of groundwater pollution, 
extent of groundwater pollution and remedial measures to be taken 
in order to rectify the effects of such pollution. 

The major sources of groundwater contamination are given 
in Table 1. 

After identifying the sources of groundwater pollution, 
the extent of such pollution must be determined. To predict 
this extent the mass transport of pollutant in groundwater must 
be studied. Several mathematical models have been described, 
that simulate pollutant; transport in aquifer. The phenomena of 
pollutant transport have been described in the following sections. 

3.1.1. Mechanical Dispersion : 

As the flow takes place, pollutant gradually spreads and 
occupies larger portion of flowr domain. This spreading is in 
both longitudinal direction, namely that of the average flow, 
and in the direction transverse to the average flow. This spre- 
ading caused by velocity variations at the microscopic level is 
known as mechanical dispersion or convective dif fusion. 



TABLE 1. CLASSIFICATION OF SOURCES OF GROUNDWATER CONTAMINATION (Lindorff, 1979) 
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3 . 1. 2. Molecular Diffusion : 

An additional mass transport phenomenon occurs simulta- 
neously with mechanical dispersion. This occurs due to variations 
in pollutant concentration within the liquid phase. This is caused 
by molecular diffusion. This molecular diffusion is of two types. 

3. 1. 2. 1. Ordinary Molecular Diffusion : When a fluid is consi- 

dered as continuum, the velocity of small fluid element is defined 
as the ratio of total momentum (vector) of its molecules to their 
mass. Actually, in addition to this velocity, the molecules in 
this element have random motion. For example, a fluid element is 
said to be 'at rest' if its centre of mass remains stationary 
(i.e. if the total momentum of its molecules is zero), although 
the molecules are moving randomly in various directions. As a 
result of the random motion, molecules of certain material in 
hi^i concentration at one point will spread with time even in 
fluid 'at rest'. This net molecular motion from a point of higher 
concentration to one of lower concentration is called ordinary 
molecular diffusion. 

3.1. 2. 2. Turbulent Diffusion : Fluid flows in nature are usually 

turbulent. The velocity at a point in such a flow fluctuates 
randomly in direction and in magnitude. A turbulent flow is often 
imagined to be a main flow with numerous eddies of various size. 

When a turbid ent flow is described with the temporal mean 
velocity q and the temporal mean concentration c at each point, 
the diffusion equation for dilute solution is 



11 


3c 

at 


+ q 




D 


m 


c 


C + V 


(Dj. V c) 


(3. 1) 


where D t = coefficient of turbulent diffusion or eddy diffusion 

and D s= coefficient of molecular diffusion, 

m 

Usually D^. >> and the term D m can be neglected. In a 
channel, the turbulence is anisotropic, especially near the solid 
boundary. For anisotropic turbulence, is not a scalar quan- 
tity and is often assumed to be a second rank tensor, with nine 
component 
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(3. 2) 


The components D- . vary from point to point in a flow. For an 

isotropic turbulent flow, = D zz and other D ± . are zero 

M D — D as D as D = D = D v „ = 0) . 

'- L * e * ^xy yx yz zy xz xz 

3.1.3. Hvdrodvnamic Dispersion : 

This term denotes the spreading (at the macroscopic level) 
resulting from both mechanical dispersion and molecular diffusion. 
Molecular diffusion alone does take place also in absence of 
motion . Because molecular diffusion depends upon time its effect 
on overall dispersion will be more significant at low velocities. 
It is molecular diffusion which makes the phenomenon of hydro- 
dynamic dispersion in purely laminar flow irreversible. 



12 


3.2. Basic Differential Equations : 


The basic differential equation describing the process 
of transport of pollutant in a saturated porous medium, and accoun- 
ting for the effects of convection, hydrodynamic dispersion, 
adsorption and radioactive decay can be expressed in its most 
general form as (derivation is in Appendix A) 


n ff + !l - n) if - iz l> D ij ft: - u i°] - d nc + (1 - n)s l + « c * 

J (3.3) 

in which 

c = concentration of dispersing mass in liquid phase i.e. 

the mass of solute per unit volume of solution (ML 
s *= concentration of dispersing mass in solid phase i.e. 

the mass of solute adsorbed per unit volume of porous 
media (ML j 

= porosity of the porous media i.e. pore volume per unit 
volume of porous media (L°) 


n 


ij 


th 


u. 


2 — * 1 

s= hydrodynamic dispersion tensor (L T ) 

= the component of Darcy* s velocity vector in the i 
direction (LT 1 ) 

= radioactive decay coefficient equal to reciprocal of 
mean life time of radioactive solute tracer (T 1 ) 

= volumetric fluid injection rate at the source per unit 
volume of medium (T ) 

= concentration of source fluid 
= time ( T) 

30 ^, Xj = the spatial co-ordinates (L) 


c* 

t 



13 


Assuming a linear equilibrium adsorption, concentration 
of solid phase can be expressed as 


s = kc 


(3.4) 


■where k is constant. 

With this assumption, equation (3.3) takes the following form 

n fl + (l_n)k || = [ n D ij fir " u i c l - x C nc + ( l-n)kc] + qc* 

° r 

- . ( l-n)k i _. 3 .. r_ ^ 3 c ^ _ _ r. . (l-n)ki . _ 

n atL 1 + n j - ax. L n D ij a Xj " u i c J " Xnc l 1 + ~ — J + 


(3.5) 


the term (l + k ~“) in equation (3.5) is constant for given 
porous media and is known as the Retardation Factor. Thus equa- 
tion (3.5) becomes 


or 
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= 0 


(3.6) 


Equation (3.6) is the general form of equation describing 
the process of pollutant transport. 

For Cartesian co-ordinates system, the above equation can 
be written in an expanded form for the case of one dimensional, 
two dim ensional and three dimensional dispersions as shown below. 

. ^ One dimensional dispersion in one dimensional flow field 

(U *°> - o 

(3.7) 
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(ii) 


Two dimensional dispersion in two dimensional field 
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(iii) Three dimensional dispersion in three dimensional flow field 
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While writing the equations in the above forms, the components of 
dispersion tensor are treared constant, an assumption which is 
valid for homogeneous and isotropic porous media. The components 
of dispersion tensor for isotropic porous media can be worked 
out in terms of dispersivity and Darcy velocity components, as 
per follovdng equations (Bear, 1972) 
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in which and D T are given by 



(3. 15) 


(3. 16) 
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(3. 17) 


where cc^ = longitudinal dispersivity 
Orj, = transverse dispersivity 

O __ -1 

= longitudinal dispersion coefficient (L T ) 

2 _ i 

D t .= transverse dispersion coefficient (L T ) 

3.3. Finite Element Formulation : 

3.3. 1. One Dimensional Dispersion in One Dimensional Flow Field : 
The basic differential equation describing the process 

of transport in saturated porous medium and accounting the effects 

of convection, hydrodynamic dispersion, adsorption, and radioactive 

decay can be expressed in Cartesian co-ordinare system 


R — D 9 — 4, — 2. — ( -n c ) + AJR C — •£ C* ss 0 
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dx 


(3.7) 


Initial and boundary conditions are: 
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(3. 19) 
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(3.20) 


Now we approximate the unknown concentration c and also the compo- 
nents of spatially dependent parameters in terms of basic func- 
tions Nj (x, y) by 
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(x, t) 


n 

r n . (x) 

j-i 3 


V. (x) 


= Nj • • • J 


9^ t) 

<Mt). 

<P 3 (t) 


t V t5 J 


= L»j J f'PjU) I 


(3. 21) 


(3. 22) 


The approximating integral equations are obtained from Galerkin- s 
scheme by making the residual, e te) . generated by introducing 
(3.21) into (3.7) orthogonal to each of the N basis function N 


/ dx = 0 


J 

(3. 23) 


(e) 


= R 


dc 


(e) 


at 


D jVf l + 1 A. (u c (e>) + \Rc te) 
D xx ax 2 n 3X x 




(3. 24) 


thus equation (3.23) becomes 
h <3 _(e) 


r. 3c - — — 

: N t L R it D *x ,2 


o 

i#e 

h 

J 
o 


3 2 c i + 1 3_ ( U c (e) ) + \Rc te) - f c^dx 

2 n 3 x x n 

3x 


(e) 


_ h aV el ^ + l K . 2- ( U c Ce) ) dx 

/ q R ft— - D xx / 3X 2 nt 1 9X x 


+ j- h N. KRc (e ' dx - / »i f ° ^ 0 

o 1 ° 


x*e. 
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f n. r dx - d n . 

j i 


3t 


xx 1 ax 


h h _(e) 

+ r D N ! 

J VY 


o o 


2 o. 

xx i 3 x 


dx 


.(e) 


h 


(e) 


+ n / N i (u x C } d * + / N i XRc " / N i n 


2. r* 


c dx = 0 


i.e. 


39 


ne 


/ [R N ± N ] dx {-^3 


D N. 




(e) 


xx i 3x 


h h 

+ Td f ( N ! N '. ) dx 
o L xx o ; 1 J 


u h .. a „ 

+ ~ / (N.N') dx + XR / [ (N.N . ) dx] {9 3 - / N ± £ c dx == 0 

n O J- J Q 1 J O 

(3. 25) 


h 


writing equation (3.25) in matrix form 
h i s h 

/ R {N3IN ! dx £93 me; + [D / {N ’ 3[N 1 J dx + 
o ° 

u h h ( np \ 

+ / {N 3 [N 1 J dx + XR ; {N}|Nj dx {9 V ' 


_ D n. -as 

“ u xx 1 i 3x 


h h 

+ -^ c* / { N 3 dx 

U * 

O o 


[p] t lf i lne) + [K] (na> !<p} (ne) = !F! (ne) 

3r 


(3.26) 


where 


h 


fp] = / R { N 3 [ N J dx 


(3. 27) 


[K] = 


h u h 

D / { N * 3[ N* J dx + — ^ / £ N 3 [ N 1 J dx 

XX 7 n r, 


h 

+ XR f £N 3 L Nj dx 
o 


(3. 28) 


£F 3 = 


D N . |§ 
xx l 3x 


h h 

+ -2 c* f £ N 3 dx 
n 


(3. 29) 
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In order to solve this equation in time domain, an impli- 
cit finite difference scheme is used for approximating the time 
derivative as given below 

{ lt } = ft ( tCp}t+At - {<p}t ) (3.30) 

For better accuracy of results column vector {9}, in the term 
[k] {9 } on left hand side of equation (3.26), can be expressed 
as an average value at time steps t and t+At i.e. 

{<p} = | ({<P} t+At + {Vf) (3.31) 

Now from (3.26), (3.30) and (3.31) 

C l_ [p] + I[ K ]) = (~[P]-f[K]) + (3.32) 

Equation (3.32) thus represents the general assembly of set of 
F.E.M. equations. Solution of equation (3.32) after proper incor- 
poration of initial and boundary conditions, will yield the unknown , 

values of concentrations at various nodal points. 

. • ) 

According to Zienkiewicz (1979), equation (3.32) can be j 

written as 

( i. [P] + epc]) = (^[p] - H-e) [k]) ITJ* * tF! j 

(3.33) | 

Equation (3.32) is the special case of equation (3.33) J 

for 6=4 (Crank-Nicholson case). J 

. ^ ■ ■ . 

Case l: To find elemental matrices for linear interpolation poly- 
~ nomials 

Let the domain be subdivided into a number of parts 
called finite elements. Let a typical element have nodes J and K. 
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Assume the solution over the element as 

c (e) * a + bx (3.34) 


where a and b are functions of time. 


At x 

= 0 

(i.e. at node j) , 

let 

c (e) = c P 1 (t) 

At x 

= h 

(i.e. at node k ) , 

let 

c (e) = 92(b) 

Then 

f 2 

= <P 1 + bh 

b = 

<P 2 - 9, 
h 


Substituting the value of a and b in equation (3.34) and keeping 
in mind that 9« and 9 2 are the functions of time. 





x 





= |/N(x) J *9(t) 3 (ne) 


where N(x) = L^tx), N 2 (x) J 


J9(t> i (ne) 


f 9 1 ( t) ] 

(i 

1 1 ~ 
[9 2 (t)j 

N^tx) = 

1 - 

X 

h 

N 2 (x) = 

X 

h 



(e) 

For linear el anient compatibility of c is needed, and complete- 
ness of c and c 1 are also needed* In our element, compatibility 
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and completeness requirements are satisfied. Now 


[p] (e) = /“ R {NJ 'LNJ dx = ^ 


2 1 
1 2 


[ K ] 


h u h 

(e) = D {N , }LN , jdx+~ / {NltN'Jdx 

XX ^ 


+ XR / InH N J dx 
o 


D 


XX 

h 


1 -1~ 

u 

x 

-i i 

A XRh 

2 1 

-1 1_ 

+ 2n 

-i i _ 

+ 6 ~" 

1 2 


{F } (ne) = {D^ N ± f§ 


(e) 


} + f c* / i N 5 


dx 



SL c * 

2n 


hi 

lx) 


. ^ e) 
f r 3 c 

1 

D xx 3 x 

J 

\ . (e) 


1 D xx 3x 

k. 


C * 

2n 


1 

1 


For one element equation (3.26) becomes 

at 


Rh 

T ! 


2 1 
Ll 2J 


39 . 

3t 


D 

“ 1 

-1 

u 

X 

"-i 

i 

, XX 

+ 

! 

- 1 

' 

3J 

+ *=- 
2n 

t > 

1 

i _ 


+ XRh 
6 


2 1 
1 2 J 



'-D ^ 


(e) 


xx 3x 


3c 


Ce) 


xx 3x 
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For twD elements. 




( 3 . 36 ) 
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For three elements 

Assembled equation will be as follows: 



Initial and boundary conditions require to be incorporated at 


3c (e) 

first node -g^r - = 0 and g x 


= concentration gradient 


At last node 


= 0 , 


and 9, 


= 0 


Thus equation (3.37) becomes 
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From equation (3.38) we will get two simultaneous equations in 
and ^ 3 . We can solve for 9 2 and C P 3 . 

Case 2: Elemental matrices for quadratic elements 


To find elemental matrices for quadratic elements, quad- 
ratic interpolation model for solution over element can be expressed 
as 

c (e) (x, t) = a 1 + a 2 x + a 3 x 2 (3.39) 

where a., a_ and a_ are the functions of time. Shape functions 

JL JL >5 

(N) are given in Appendix B. 

[p]^ = R / l N 2 L N j dx 

o 
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Again 

, . h u 

[K ] teJ «= / [d {N , UN , J+~{NHN'J-tAR{N}[Nj]dx 

“ " ^ xx n 

o 
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3 9 . 

Now at first node 5 - — 

o r 


26 



Thus equation (3.41) becomes after incorporation of initial and 
boundary conditions 



Equation (3.4 2) can be solved for 9 2 # 9 3 and ? 4 



27 


3*3.2. Two Dimensional Dispersion in Two Dimensional Flow Field ! 
Governing differential equation 


3c r « Jc , „ 3 2 c . „ 3 2 c . ^ 3*c 


R H-[ E xx^1 +D xy 3x3y 


+ D 


yx 3y 3x 


+ D 


yy ~~2 
3y 


1 


+ n th (u x c) + fy (u y c)] + XRc - f c * = 0 


Initial and boundary conditions are: 


(3.8) 


y, t) 

as 

c o 

for 

all 

X 

«* 

*i 

1 V 

0 

and 

t *= 0 

(.3. 43) 

y, t) 

s 

c o 

for 

t > 

0 and 0 

< 

Y < 

0 

(3.44) 

y* t) 

s 

0 

for 

t > 

0 and 0 

< 

Y _5 

0 

(3.45) 

x, 0, 

t) 

= 0 

for 

t > 

0 




(3. 46) 

x, 0, 

t) 

= a 

for 

t > 

0 and x 

> 

0 


(3.47) 


3y 

_3c 

3y 

Let the trial solution of equation (3.8) be assumed of the form 

n 

j 

K (t >] 


c(x, y, t) 


c (e) (x, y, t) = S [K.(x, y) J C«P(t) ? (ne) 

i =1 J 


= lN r N 2 ... N n ] 


Cp 2 ( t) 


9 ( t) j 

n J 

(3.48) 

where EL (x, y) are the interpolation functions. <P(t) are the 
unknown nodal values of concentrations. 

Ihe approximate function c^ e) (x, y, t) will be exact solu- 
tion to the equation (3.8) if the error or residual 
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(e) 


= R 


3c 


Ce) 


3t 


r a 2 c (e) ( e ) *2 (e) a 2 (e) 

■ f D xx ^2 + D xy 3 y 3x + D yx 5^ + D - 


_yy 37" 


+ n (u x c(e)) + fy ( V (e))1 + XRc(S) ' n c * 
is equal to zero. 

To minimize the residual, applying Galerkin method, 

(e) 

/ / N. e ; dA = 0 

A (e) 


(3.49) 


11 N i ft" 

x(e)' 1 3t 


(e) 


22 (e) -2 (e) - 2 (e) a 2 (e) 

- (D ... M + D ... 4^ + D ... + D -2-S ) 


xx 2 
3x 


xy 3 x 3y yx 3y3x yy a 2 


+ n^hi (u x c(e)) + % (u y c(e) )l + ^Rc (e) - g c*]' dx dy = 0 

(3.50) 

-2 (e) r 3N.. a _(e) 


r f D N. *“-5 
* , \ xx 1 . - 

A (e) 3 x 


*2 (e) 

/ / N, D r 

A (e) 1 ^ 3y 2 


N. 

3c 

1 

ds 

- / 

/ 

l 

3x 

X 

A 

(e) 

N. 

3c 

1 

ds 

- / 

f 

1 

3x 

X 

A 

(e) 

D 

N . 

3c 

1 

ds 


yy 

1 

3y 

y 



3N. 3N, , v 

5 l. _J. ^{cp } (ne ' 
xx 3x 3x 


3N. 3N, , v 

/ / D dA {®} (ne) 

A (e) yy 3 y3y 


3 2 c 


r r D N u — - dx dv = (£ D N . ~ 1 ds 

J\ { xy 1 3 x 3 y y J xy 1 3 y x 

A e 


3N. 3N. 


- / / D 


(e) 


xy 3 x 3y 


dA {<P 3 


(ne) 


/ / D N. 

(e) yX 1 3y9x 


9 " ■ - dx dy 


s= N.. — 7 : 1 .. ds 


3e 

'yx ' i ax ~y 


3N. 9N -i 

- / / D —A — j- dA if i ne) 

A (e) yx * 3x 
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/ f — N_, ~ dx dy = 


(e) 


n i 3x 


r K i cl x ds -' ( ' ^ i<p!<ne> 


U 

/ / N . - 5 ^ dx dy = 

(e) n 1 ^ 

A Ve; 


u 

— ^ N . cl ds j j 
n x u x y (e ) . n dy 3 

£\ 


u 3N - / __ ^ 
/ / — -± N. dA {9 } ( e) 


Now equation (3.50) can be written as 


A 


(e) 


1 J 


99 

3 1 


(ne) 


3 N . 3N. 
1 1 


3 N . 3N . 
1 . _1 


A 


(e) 


3N, 3n. 

4 . q ■ J. j. 

yx ay ax 


D 


xx ax 

ax 

<r jl> ~ 

xy 3x 

3Kb 

a Nj 

u 3N * 

x i n . 

yy ay - 

ay 

n ax j 


u 3N . 

__Z — £ n . + AR N . N . ] dA i 9 } 
n ay j i J 


(ne) 




( D IS D l£ _ Is c ) i + (d — + D — - — ^ c) 1 IN. ds 
+ u ' “x v vx a x vv av n y-J i 


xx ax xy ay n 


+ / / - c*{ N . 3dA 
/ \ n x 

A e 


(3.51) 


Writing equation (3.51) in matrix rorm 

/ / R { N 1{ N J dA {9 } (ne) + / / (D xx tN, x 5 [N /x J + D xy {N, x }[.N, y ] 

A (e) A (e) 

+ D yx {K,y)lN, x 1 + B yy fN, y }|N, y J - iN <x }[N J - ^ f%3L N ] 

+ \R IN }{ N J) dA (9 S (ne) 

ir/r* 3c xn is-ds c )l + (D — + D - - ^ c)l 1 N. ds 

= jL (D xx3lc D xy 3y n c ;i x 1 yx ax yy a y n y J i 

(3.52) 


+ / / c* { N } dA 

A (e) n 


y 
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Equation (3.52) can be written as 

[p] {q>} (ne > + [ K ] {<p} (ne > = { F B } + { F } 


(3.53) 


where 

[ p ] 

[k] 


= / / R CnJLn J dA 


A 


(e) 


\i) <Dxx fN 'x !lN 'x J + + D vx J 


xy 


yx 


{F b } - 


+ D yy W.yUM.yJ- 5» In, x !Ln J - JE !N, y 3lnJ 
+ AR {N }[n]) dA 

§ [ <D xx ll + D xy % - =* o)1 - 


xx ax xy 


+ (D ~ + D — 
yx 3x yy 3 y 


u v 

- x ) i In. ds 
n y J i 


{F} = / / ^ c { N } dA 


(e) 


n 


In order to solve the equation in time domain, an implicit finite 
difference scheme is used for approximating the time derivative 
as given below 


{<P } 


(ne) _ {<P3 t+At - j? f 

At 


For better accuracy of results, the column vector { $ }^ ne ^ in term 
[k] {9 on left hand side of equation (3.53) can be expressed 

as an average values at time t and (t+ t), i*e* 

{<p} (ne) = \ ( {<P 3 t+At + 19 3 t ) 

Substituting the values of {<p}^ ne ^ and }^ ne ^ in equation (3.53) 



31 


(•— [p] + \ C K]) C«P l t+At = (A- [p] _ 1 [K]) £9 f + {F 3 + {Fg ] 

(3.54) 

Equation (3.54) thus represents the general assembly of set of 
F.E.M. equations# solution of equation (3.54) after proper incorpo- 
ration of initial and boundary conditions, will yield the unknown 
values of concentrations at various nodal points. 

Now according to Zienkiewicz (1979), equation (3.54) can 
be written as 

(~ [p] + e [k]) {<P } t+At = [p] - (l-e) [k]) {<P f- + {Fg 3 + £F 3 

(3.55) 

Equation (3.55) is the special case for © = -| (Crank-Nicholson 

case) . In order to solve the equation we need compatibility of c 

and completeness of c# and For this assume a suitable form 

gx gy 

(e) 

of variation of concentration c inside a finite element 1 e' . 

We assume a linear variation as 

c te ^ (x# y# t) = a 1 + a 2 x + a 3 y (3.56) 

where a 1 # a 2 and a 3 are the functions of time. 


Calculation of Elemental Matrices: Shape functions (N) for tri- 


angular elements are given in Appendix B* Now 


3tl l 

b l 

3N2 _ 

b 2 

3H 3 b 3 

3N 

3 x 


dx 


3x " 2A^ e ^ # 

w 


3N 2 c 2 

3N 3 c 3 

3Y 2A (e) ' 

3y “ 
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Now 

[p] 


/ / R { N }[ N 1 dA 
*(e) 


/ / R 
»(e) 


N. 

n] 

« 

N. 


\ 


IN. 


N, 


N-] dA 


/ f 

Ae) 


R / / 

A (el 


RA 


(e) 


12 


» 2 1 

N„N_ 

1 2 

N 1 N 3 


N 2 N 1 

N 2 

N 2 N 3 

. dA 

n 3 Ni 

N 3 N 2 

N 3 2 


1 

1*1 


■ 


2 

1 


D xx' N -x ! L k -xJ ® 


1 

2 


D / / 
xx , x 

A Ce) 


■f 


3N ^ 

3x 

3^ 

3x 

3N_ 


l 3x 


3 N. 


3N„ 3n. 


3 x 3 x 3x 


J dA 


D / / 


2A 


b l 

liT 


2A 


l 2A 


(e) 


V 


Ce) 


2A 


eK 


b 3 J dA 


D 


4A 


xx 

(e) 


"1 

b ^l 

b 3 b l 


b l b 2 


b 3 b 2 


b l b 3 

b 2 b 3 
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c 

c 

c 

/ / XR CN 3 LN.J dA - 12 

A le) 


/ ' ^ tN.yJlNJ ^ = 6? 

A Ce) 
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D 2s + D |s - 2s c)1 _ 

D xx 3x xy 3y n 


(D + D 

V -TTV A 'V \ 


U 

3c _ _y 


V , 1 


on boundary vfiere §“ 





tH CM CM 
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r 





. .X 




) + a D c)l 

xx 9 x xy n 


+ (D vx fS + 


n n 


u 

v^ , 


Let the 1-2 is boundary of an element triangle 12 3 . 
Let distribution on 1-2 is 


N. 


N, 


N cp + n cp 

ri 2 y 2 

N. 

x h ^ s 


1 

h 


N. 


1 

h 


dy 

dx 


= tan©. 


dX = 
ds 

dx _ 
ds ” 


sin© 

cos© 


- 3 N i 9 N o 

2 c _ ± cp + £ Cp 

9 x 9 x 1 3 x 2 


_ 3 s cp + ^2 9 s ^ 
3 s 3 x 1 3 s 3 x 2 


_as 

3 y 




JL Cp -H — cp 

hcos© T 1 hcos© T 2 


3N 9N 

i <p 4 . _ — £ cp 

3y Y 1 5y 2 


hcos© ^ 9 2 " ^1 ^ 
9 N. 


9 N 9 N , 

1 3-§ cp 4 . — m _ -i 

3s ay 1 3s ay 2 hsin© 


^ 2 ^ 1 ] 


D xx fx N i ds = q D xx hcos 0 (Cp 2 _Cp l ) 


s l 
_ h 

s 

h J 


i x ds 
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n 




' Cheese 'W 


fl_ J 1 
h 

D ! 

'h/2 

I s 

1 h 

r cose ds - (<p 2 -<p 1 ) ( 

[h/2 . 


D 


(m _cp ) 

2 v z Y r 


I 


1 1 
1 

h 


D 


xx 


D cl N. 1 ds = / D a cos© 


xy xx 


xy 


D 

= — — ^ a hcosS 


-1 1 

-1 1 

s 


5 

h J 



ds 


Kx^ 1 y N i dS ' C / D yx hcosS (tp 2- ,, l ) < 

h/2 


1 - 
1 h 


s 

h 


sine ds 


= tan© -jpS j 


— tan© 


r-i i 
-1 1 


h/2 

(9 


<P. 


fas D 1 N . ds 

J yy ay y x 


h 

/ D a sin© 

yy 
o J - 1 


c N. i x ds = / hs LVl. + 


1 - -S' 1 

1 h 


_s 

v. h 


D 


V ds = — 5p a sin© h 



N. 


l N 2 


cos© ds 


u h cos© 
X 

6n 


h u cos© 

-Jv 


on 


+ *2 

*1 + 2(P 2 


2 1 
1 2 


«P. 

9 


2J 
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Similarly 


$ % M 1 . h U Y sinG 
J n c ds = 


6n 


‘2 1 " 


I’M 


1 


_1 2 _ 

1 

L 9 2J 


[' 8 ll - 


. |S + D a - ^2 C )1 

xx 9x xy n x 


+ (D, 


3c 


u 


+ a " 7? C) Kr 1 N -! dS 


yx 3x yy “ n w "y J "i 


= (■ 


xx 


1 1 
1 1 


+ tan© 


- 1 1' 
-1 1 


h u x cos© 

_ — 


2 

1 


r 

2 




sin© 


6n 


2 1 
1 2 


t m •) 


2 J 


+ -rp' a h cos© 


.1 


D 1’ 

+ — o. h sin© 

2 , 1 . 


= [ H§* + ^ tane) 


■ 7 — (u cos© + u sin©) 
bn x y 


- 1 1 
-1 1 
"2 1 
.1 2 J 


<?. 

<P. 


+ | h a (D^ cos© + sin©) ■) 





® *►) 
c _>-> ) 


h = 


((Xj-X ^ 2 + (y 2 -y 1 ) 2 ) 1/2 


(b 2 - -V/2 
vo 3 ' w 3 / 


tan© = 


sin© 


cos© = 


?2-r i 

x 2 -x l 


_^3 _ 

2^1/2 - 


(b^+cj) 
c 3 

(b|+c|) 1 / 2 
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M - t-2 (D 


-(d - d — ) 

2 v XX yx c 3 ; 


-1 IT 

L-i iJ 


+ (u c 0 - u b_) 
on x 3 y 3 


2 1 
1 2 


+ 2 (D xy c 3 D xy b 3^ 


Similarly when 2-3 are the boundary nodes 




D D b 

xx yx _lv 

2 2 c 1 ; 


-1 1 

-1 1 


% (d. 


c - D b „ ) 
xy 1 yy 1 


+ -r~ (u c -u b.) 
6n x 1 y 1 


f2 1 1 

] 

JM 

- 1 2 J 

9, 

l O J 


When 3-1 are the boundary nodes 

1 


[ f b) 3 - t| (D : 


XX 


b 2 
yx c ' 


-1 


1 

1 


+ 6H (u x C 2- u y b 2 } 


'2 1 
1 2 


] 


9, 

9. 


+ 2 ^ D xy c 2 " D yy b 2^ 


Equation (3.8) is also solved on the same lines as above 
for second set of initial and boundary conditions which are given 
below: 


c(x. 

Y. 

t) » 

c o 

for 

all 

x, y > 

0 

at t = 0 

(3.43) 

c (0, 

y# 

t) = 

0 

for 

t > 

0 and 0 


Y < P 

(3.57) 

c( , 

Y> 

t) = 

0 

for 

t > 

0 and 0 

< 

Y < P 

(3.45) 

ff lx - 

o. 

t) = 

0 

for 

t > 

0 



(3.44) 

c (x, 

Pt 

t) = 

0 

for 

t > 

0 



(3.58) 
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CHAPTER 4 

RESULTS AND DISCUSSION 

4. 1. General : 

In the present work one dimensional and two dimensional 
dispersion equations for predicting the pollutant concentration 
in groundwater flow have been solved using the finite element 
technique. In the one dimensional case both linear and quadratic 
elements have been considered. In two dimensional case triangular 
elements- have been taken and the equation has been solved for two 
different sets of boundary conditions. In solving a dispersion 
equation a suitable value of © has to be cnoosen. To find the 
optimum 6, the one dimensional equation was solved for different 
values of 0 and time intervals. It was noticed that at the value 
of © = 0.75 the analytical results matched closely with the 
results obtained by the finite element method. The concentration 
values for different values of 9 ror t = 15 min, are shown in 
Table 2. Based on these results © is cnoosen to be 0.75 for all 
subsecruent calculations. 'The various results for both one dimen- 
sional and two dimensional cases have been discussed in the 
following paragraphs. 

4. 2 . One Dimensional Case. 

t, 

In one dimensional case dispersion equation 

was solved by taking both linear and quadratic elements and results 
obtained are compared with the analytical results. 
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4 . 2 . 1 . Linear Elements : 

Case 1 * The dispersion equation was solved for the 
special case where the retardation R is taken equal to 1. Volu- 
metric fluid injection rate q and concentration of source fluid 

* 

c are respectively zero velocity in x direction is taken as 
0.039 m/min and porosity as 0.39. The equation was solved by 
using both analytical and F.E.M. methods for different time 
intervals. The results have been presented in Figure 4.1. It is 
seen that the results obtained by the two methods match closely 
with deviation never being more than 3 percent. Slight deviation 
in the finite element results might be due to the following 
reasons 1 firstly inaccuracy may result due to semi-implicit- 
explicit finite difference scheme used for approximating the rime 
derivative. Secondly, due to comparatively larger value of time 
interval used. Thirdly, due to the effect of finite boundary 
which was taken as 50 m in this case, lastly, due to the length 
of element choosen. 

Case 2 ‘ The same equation was solved by changing the 
value of R from 1 to 1.156. The corresponding results are shown 
in Figure 4. 2. 

Case 3 ‘ In this case the value of X (for radium) is 

changed from 0 to 4.343 x 10 _S per min. Tne results are found to 
be the same as for Case 2. Thus it may be concluded that when 
the value of X is less it does not have significant effect on the 
concentration distribution. 
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Case 4: 


In tnis case trie value of q is changed to 0. 006 


_ , * 

per mm ana c to 10 rag/L and all o ther values are same as in 
rirst case. The results are shown in Figure 4.3. .-it some dis- 
tance rrorn the source, the concentration is more than initial 


concentration. This is due to volumetric fluid injection rate 
and concentration or source fluid. It is seen that after some 


distance concentration stabilises. 


Case 5: In this case, the value of R (retardation factor) 

is changed rrom 1 to 1.156 and all other values are kept same as 
in previous case. The results are shown, in Figure 4.4. It is 
seen that the nature of variation is same as in Figure 4.3, 
though because of adsorption pollutant concentrations are less in 
this case. 


Case 6: In this case all values are kept same as in 

—8 

first one. Only the value of a was taken as 4.34 x 10 per min. 
This variation is shown in Figure 4.5. 


Quadratic case : The dispersion equation was solved for 

all the cases described in previous sections using quadratic ele- 
ments. The results match closely with those obtained for linear 
cases (Figures 4.1 to 4.5). 

Comparison of results when the problem was solved by 
linear approximation (Figure 4. 1) and quadratic approximation 
(Figure 4.6) for the same set of values show 1 that the results 
closely match each other. It is also seen that the results obta- 
ined by quadratic method match even more closely with the analy- 
tical re.sul ts than those obtained oy linear method* 




42 


4.3. Two Dimensional Cssp : 

Two dimensional dispersion equation is solved using finite 
element method by using triangular elements. The governing equa- 
tion was solved for two different sets of boundary conditions 
described below. Results for one set of boundary conditions are 
given in Figures 4.7a and 4.7b. 

The domain was taken as rectangular area 100 m in x direc- 
tion and 30 m in y direction.. This area is divided into square 
elements 5 m x 5 m which are again divided into triangular elements 
thus giving 240 elements with 147 nodes. Initial concentration of 
pollutant is taken as 5.0 mg/L all over the domain. In this case 

the value of retardation factor (R) is taken as 1.0, D and D 

xx yy 

*“5 2 “~6 2 

are taken as 6 x 10 m/min and 9 x 10 m /min respectively, 
the velocities in x direction and y direction (u x and u y ) are 
taken as 0.06 m/min and 0.012 m/min respectively. Value of cl is 
0. 1 mg/L per meter and the values of cross coefficients of dis- 
persion D yx and X, q and c are all zero. 

The dispersion equation is solved for the first set of 
boundary conditions shown in Figure 4.8 The concentration at 

(different time intervals are found at all the nodes in the domain. 
From this, contour of concentration of pollutants at various time 
intervals are plotted. 

Figure 4.7(a) shows the results for t = 10 min and Figure 
4.7(b) for t = 25 min. From Figures 4.7(a) and 4.7(b) it is seen 
that near impervious boundary concentration gradient is high. 

Near the boundary at vhich the concentration is specified as 



5.0 mg/L, the concentration gradient is also hi^i. As we go 
away from the impervious boundary as well as specified boundaries 
the concentration is almost same (i.e. concentration gradient is 
very small) and it decreases with time. 


TABLE 2 . COMPARISON OF RESULTS FOR LINEAR ELEMENTS AFTER 5 MINUTES 
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Table 2 (continued) . COMPARISON OF RESULTS FOR LINEAR ELEMENTS AFTER 10 MINUTES 
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Table 2 (continued) . COMPARISON OF RESULTS FOR LINEAR ELEMENTS AFTER 15 MINUTES 










Table 2 (continued). COMPARISON OF RESULTS FOR QUADRATIC ELEMENTS AFTER 5 MINUTES 
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Table 2 (continued). COMPARISON OF RESULTS FOR QUADRATIC ELEMENTS AFTER 10 MINUTES 
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Table 2 (continued). COMPARISON OF RESULTS FOR QUADRATIC ELEMENTS AFTER 15 MINUTES 
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CHAPTER 5 

SUMMARY, CONCLUSIONS AND SUGGESTIONS FOR FURTHER WORK 

5.1. Summary : 

In water resources development and management one of the 
main problem is that of water quality. While discussing the water 
quality problem the question of transport and accumulation of 
pollutant in an aquifer needs special attention. Though it seems 
that groundwater is more protected against pollution than surface 
water but when groundwater is polluted, it is very difficult to 
bring it back to its original state. 

Most of the groundwater becomes polluted by the pollu- 
tants brought down from the surface. Under natural conditions an 
equilibrium exists between the water leaving the formation and 
pollutant which is carried with this water. When an aquifer is 
recharged with an inferior quality water or with water having 
more pollutant concentration than that already present in the 
aquifer, this equilibrium is. destroyed. The equilibrium is 
restored again after a rise in pollutant concentration in the 
aquifer. Sometimes this new pollutant concentration may go beyond 
the permissible limits, thus causing pollution problem. 

The sources of groundwater pollution may be the following, 
of water through carbonate rocks, sea water intrusion into 
aquifer, accidental breaking of sewers loading to discharge of 
inferior quality water into the aquifer and percolation from 
septic tank. Sometimes rain water and irrigation water carrying 
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the fertilizer salts, herbicides and pesticides may infiltrate 
through the ground surface polluting the aquifer. After studying 
the various causes of pollution of groundwater it becomes essen- 
tial to have methods for prediction of the distribution of the 


pollutant in the groundwater. 

In the present study an attempt has been made to study 
the pollutant transport in groundwater flowing through porous 
media. The spread of pollutant in groundwater takes place in 
many different ways like mechanical dispersion, molecular diffu- 
sion, hydrodynamic dispersion etc. Also other phenomena like 
adsorption, radioactive decay etc. may affect the pollutant trans- 
port phenomena. So suitable mathematical models need to be 
developed including all the above mechanisms in order to suitably 
represent the actual conditions. Keeping this in view a suitable 
mathematical model has been developed to adequately describe the 
pollutant transport phenomena in the groundwater. The resulting 
differential equations both for one dimensional and two dimen- 
sional cases were solved using finite element technique. 

To solve the dispersion equation for pollutant transport 

in one dimensional flow field, both linear and quadratic elements 

have been taken ror Galerkm Finite Element solution. Special 

cases of same problem were solved analytically and results were 

compared with finite element solution. The results match closely. 

Further, the solution of dispersion equation was done by varying 

* 

various parameters like R, A, q, c and results are plottea ana 
analysed. It is noted that if the value of A is low then its 
effect on concentration profile is negligible. Also the 
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volumetric fluid injection rate and concentration of source of 
fluid affect the pollutant concentration at a distance far away 
from the source. It is also seen that the results obtained by 
linear and quadratic elements give matching results. 

For solving two dimensional dispersion equation in two 
dimensional flow field, tvo different sets of boundary conditions 
were considered. Results obtained from the first set of boundary 
conditions are plotted in the form of concentration contours. 

These results indicate that near the impervious and concentration- 
specified edges, concentration gradients are high and this gradient 
increases with increase in time interval. In the middle of the 
domain the concentration is nearly constant. Its value decreases 
with increase in time. 

5.2. Conclusions * 

The following conclusions may be drawn from the results 
obtained by solving one dimensional and two dimensional pollutant 
transport equations by finite element method. 

(i) During finite element formulation of the problem, the 
value of © was varied over a range of 0 to 1. It is 
found that for © = 0.75 the finite element results match 
closely with the analytical results. 

(ii) For radioactive pollutant of higher half life period like 
radium etc., radioactive decay does not affect signifi- 
cantly the pollutant distribution in the aquifer. 
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(iii) An increase in retardation factors (considering adsorption) 
leads to lesser over all pollutant concentration in the 
aquifer. 

(iv) Groundwater recharge with pollutant concentration will 
have considerable effect on the over all pollutant con- 
centration in the aquifer. 

(v) In two dimensional case the concentration contours vary 
greatly depending on the boundary conditions with concen- 
tration gradients vary high near impervious and concent- 
ration specified boundaries. In the middle portion of 
the domain in two dimensional cases, the concentration 
has been found to be nearly constant and the concentration 
gradients small. 

Some of the limitations of the study are: ( l) ion-exch- 
ange phenomena was not included, (2) three dimensional aspects of 
the study were not considered and (3) only single boundary condi- 
tions and triangular elements were used in the finiue element 
formulation in two-dimensional case. 

5.3. Suggestion for Further Work : 

Following are a few suggestions -chat can be incorporate 
in further work: 

( l) Two dimensional case results can be obtained for dif rerent 

•¥? 

boundary conditions and for different values of R, q, c and 
X. Further, 8 noded isoparametric elements can be taken to 
solve the same governing difrerential equation to give 


better results. 



Ion- exchange phenomena should be included in the governing 
deferential equations, k program can be developed for 
uhree dimensional dispersion equation. In this case in 
must be kept in mind chat this involves a lot of comparer 
memory space and execution time* 
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APPENDIX A 

DERIVA TION OF THE GOVERNING DIFFERENTIAL EQUATION 


A. 1. Pick's First Law of Binary Diffusion : 

Consider a small volume of mixture of fluids A and B. 
Let C A °B mass of A and B per unit volume of the mix- 

ture. Let and Vg be their respective average velocities. 

Then for the mixture 


c 


+ Cg and velocity V 



(A. 1) 


Thus for an isotropic material, Fick's law of binary diffusion 
relates the motion of A, relative to average motion of mixture, 
to concentration gradient 

c a (V a - V) = - c D m grad (■—) (A. 2) 

where D = coefficient of molecular dispersion of the mixture, 
m 

Similar equation can be written for material B. 


A. 2. Derivation of the Differential Equation: 




ICaUA 


I 

4 



iU * 4 * Cc*u a) 4 *. 
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Consider an elemental volume as shown above. Now 
\ = 1 U A + J v a + k 

Mass of pollutant entering the elemental volume in x direction in 
time dt = (c A u A > dy dz dt. 

Mass of pollutant leaving the elemental volume in x direction in 
time dt 

= °A U A lc A U A )dx ] dy dz dt 

Net amount of pollutant which enters elemental volume in x direc- 
tion 

= (c A u A )dy dz dt - [c A u A + ^ (c A u A )dx] dy dz dt 
= ( c a u A^ ^ dy dz dt 

Similarly net amount of pollutant which enters in elemental volume 
in y direction 

= " % ^a v A } ** dy dz dt 

Net amount of pollutant entered into elemental volume in z direc- 
tion 

= - h (c a V ^ dy dz dt 

Net increase of mass 

3 c ^ 

= at) n dx dy dz + (*H dt) ( 1-n) dx dy dz 

9t ou 

where n = porosity 

s = mass of pollutant adsorbed per unit volume of the 


porous medium. 
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Mass balance equation, when no chemical change is involved, is: 

Net mass out + Net mass out + Net mass out + Added - Mass 
flow in x flow in y flow in z mass decayed 

direction direction direction * 

= Increase of mass 

<hi ( C A U A ) + fy (c A v A } + dy dz dt + qc dx dy dz dt 

- A [ac A + (l-n)s] dx dy dz dt 

r 3c A 3 s - 

= L n g"£- + ( 1-n) -gr~ J dx dy dz dt 

where c* = concentration of source fluid 

q = volumetric fluid injection rate. 

i . e. 


- h^ (c A U A> *h '°A V A> + + < 3 C *- X t nC A + (1 - n)s l 

3 3=5 

= n W + ( 1_n) 3t 


3c^ g 

- div(c A .V A ) + qc* - A [nc A + (l-n)s] = n + (1-n) 


(A. 3) 


From Fick's law 


~ (V, - V) = - c D grad(-^) 

n A c 

Takinq D as a dispersion tensor and c average density of mixture, 
as constant 

c V = c*V - n D grad(c,) 

Now substituting in (A. 3), we get 
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- div [ c A V - n D grad(c A )] - X [nc. + (l-n)s] + 


qc 


3c 


- n rr * (1 - n > If. 


or 


ac 


n + (1-n) || = div [n D grad(c A ) - c a v] - X [ nc A + (l-n)s] 


+ qc 


Now dropping c A as c we get 


n _2£ (1 ) £s 

n at ni 3t 


div [n D grad(c) - cV ] - X [ nc + ( 1-n) s] + qc* 


In tensor notation, equation can be written as 


n ft + ll ~ n > H - air D ij -fer - cu i3 - x t nc + qc 
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APPENDIX B 

TO FIND SHAPE FUNCTIONS 

B.l. For Quadratic Element : 

Let the solution over element for the quadratic interpo- 
lation model be expressed as 

( e) 9 

c (x, t) = a 1 + a 2 x + a 3 x (B. 1) 

where a^, and are the functions of time t. Since there are 
three, constants a^, cij, a 3 in equation (B.l), the element is 
assumed to have three degree of freedom, one at each end and one 
at middle point. 


I 2 3 

K W — 

Let 

c (e) (x, t) = ^ at x = 0 

c (e) (x, t) = 9 2 at x = h /2 

c (e) (x, t) = 9 3 at x = h 

where all 9 ' s are the functions of time, t. 

From (B.l), substituting the values of 9^, 9 y 9 3 ror 

various values of x we get 
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9 3 “ ^1 + + 

a 2 = (49 2 - 39 ^ - 9 3 )/h 

a 3 = 2((p l " 2tp 2 + 


by substituting values of a^, a 2> in equation (B. 1) we get 

c (e) (x, t) = <f 1 + (49 2 - 39 1 - 9 3 )g + 2(9 1 - 29 2 + 9 3 )^ 

h 

- + C S'T5 )<P 3 

n n n 


_ (i _ 2) (i . 22 $) cp + -l2( i_ 2S) cp _ 2S( i _ 22 $) cp 
- U- h Ml - h ;H i + h^ 1 h JT 2 h ^ 1 h' T 3 


where 


= iN^tx), N 2 (x) , N 3 (x)J 


9. 


L 9 


3 J 


(B. 2) 


[ N (x) J 


{cp } (ne) 


N (x) 

Ln 1 

(x) , 

N 2 (x) , 

sz; 

M 

N 

& 

II 

(i - 

X) (i _ 2X) 
h M1 h ' 

n 2 (x) = 

i* (I _ i 
h 11 1 

5> 

N 3 (x) = 

X 

~ h 

(1 - 

2xn 

h ; 

{9 }^ ne ^ 


- cp " 
1 1 1 




1 *2 

2 i 

r ■ 



L9 3 J 



and 
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B.2. For Triangular Element: 


3 < *j,Vj 1 


The two dimensional simplex element is a straight sided 
triangle with three nodes one at each ? 
corner as shown in figure. Let the 
nodes are labelled as 1, 2 and 3 in 
anticlockwise direction. 

We assume a linear variation 
of concentration c inside ,the element 



.(e) 


c (x, y, t) = a 1 + a 2 x + a 3 y 

where CLy a 2 , a .3 are the functions of time. 
By using nodal conditions i.e. 


(e> (x r 

t) 

= c P 1 (t) 

at 

node 

1 

(x 1# y^ 

(x 2 # Y 2 , 

t) 

= V t} 

at 

node 

2 

(x 2 , y 2 ) 

(x 3 / y 3 / 

t) 

= ? 3 (t) 

at 

node 

3 

lx 3' ^3^ 


lead to system of equations 

9 1 = a l + a 2 X l + a 3 y l 

^ 2 = a l + a 2 X 2 + a 3 y 2 

<P 3 = a l + a 2 X 3 + a 3 Y 3 

The solution of above equations yield 

(a^ + a 2 <P 2 + a 3 f 3 ) 


a. 


a. 


2 A 


2A 


(Vi + * b 3 V 


(B.3) 


(B.4) 

(B.5) 

(B.6) 
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“3 

II 

if (c ih 

+ c 2 ? 

2 + c 3 Cf 

3^ 

where 

A (e) 

is the 

area of 

triangular element 

A (e) 

1 

“ 2 

(x l y 2 + 

x 2 y 3 + 

x 3 y i ~ 

y l X 2 " 

y 2 X 3 

and 







a l = 

X 2 Y 3 " 

x 3 y 2 

b l 

= y 2 - 

y 3 

C 1 

a 2 = 

x 3 y i - 

x l y 3 

b 2 

= ^3 ' 

y l 

c 2 

a 3 = 

x l y 2 - 

x 2 y 1 

b 3 

= y i " 

y 2 

C 3 


On the substitution of the values a^, a 2 , a 3 
and rearranging the terms, we get 

(x, y, t) = N^x, y) c P 1 + N 2 (x, y)9 2 + 

= |_N(x, y) J C*? C t) }^ ne ^ 


LN(x, y)J 

= [N^x, 

y) # 

V*' 

y)/ n 

N 1 (x, y) 

1 

(a l 

4_ Vn -v 

+ c x y) 

= 2A (e) 

T JJ 

N 2 (x, y) 

1 

(a 2 

4. v* 

+ c 2 y) 

2A (e> 

d 2 X 

N 3 (x, y) 

1 

(a 3 

+ b 3 x 

+ c 3 y) 



and { 9 (t) 3 


1 9 9 (t)f 

l9 3 (t)J 


123, given by 
- y 3 x 1 ) 



in equation (B.3) 


3 (x, y)9 3 (B.7) 


(x, y) J 
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c «ATN PROGRAM FOK UKEAP KlEmENTs 

c .,,„.*„ ss . ss$s;; * 5;ga; 

S2Sg£SElSSgSJ£ ,r Kf 28IH5 F3fi 

el c =elemental n*> matrix 
SLGKsSLDSAL W MATRIX 
GLFsGLOSAL (F> MATRIX 
NNOD=NlHBER OF NODFs 
nel£m=m jmber of elements 

¥?5S2 s !l!iJJJ§ Ef 5 f ,9Qr-^ on£: 5 Ax WHIC ^ concentration is spfctftfp 

W HICH C3NCF NTRA r l ON TS SPECIFIC"'* 

SPC3 = VALORS OF SPECIFIED CONCENTRATION 
UX=JARCY 5 VELOCITY VECTOR IN X-PlRgCr TU*! 

. DXXsDISPERSirm COEFICIENT in X-DTRECTlflN 
POP=PDROSITY - ^ * 

R£TA = REXAP, PATIO N FACTUR(R) 

or=riME interval 

Q~ VUL U M ETRI C FLUID INJECTION RAXES 
CSTsCDNCRNTRATION OF SOURCE FLUID 

C0MMGN/Hl/CDRD(NK0n # ?),NN3D 
C0M!K0N/A2/K0MECTlNEtEH,2) ,NELE* 

CO«MON/AB/ELGKC2 £ 2),RLF(2) 

COMMON/A4/DXX, RET A, CST,XAMn,Q,UX,POR,ur 
COMMON /A5/C0 M Cl NNjH) 

CDMM0N/A6/GLGK(NN0n,NN0 D) ,GLF(NnOd) 
COMMnN/A7/N$PCO.ISPC(2) / SpCO(2) 

COM MON /A1 /CORD (101 , 2) ,NfiOp 
C0MMOu/A2/K0NECTUO0f2)rNELEM 
COMMON/A3/FLGK(2,2) ,RLF(?) 
CQ*MD\/A4/PXX,RETA,£sT,XLA M D,0,UX,P3P,or 
COMMON/AR/CONCClOl) 

COMMON/ A6/GLGKC101 , lOl) , GLF( 1 01 ) 
COMMON/A7/NSPCO.ISPCO(2)/SPCO(2) 

DATA DXX, RETA,CST, XlAMD,0 , UX , POR , DT/O # 1 , 1 . C , 0 . ? , 0 . ° ,0 . 0 , 0 . 03« , 

1 U.39,0.5/ 

DATA NS PCO* I SPCu,SPC0/2, 1,101,10.0,0.0/ 

DATA NHJD r REL£M/.101- # IOO/ 

DO 1 Isl-fiNOD 

1 READC23,*) tC0RD(I,J),J=l,2),C0NrCT) 

DO l I=i,«EL5« 

2 READ (23,*) , (KONfiCT C 1 , J ) , 0 = 1 ,2) 

T = 0.0 

DO 3 1=1,10 
T=T*PT 
CALL MSMBL 
CALL SPCONC 

CALL GAUSJD(GLGK,GLF,NNOD) 

TYPE * , N N 0 D 
DO 49 Jal.MNnD 
CONCC J)=GLF( J) 

49 tfRITE(40,*) ,T, J r CONC(J) 

3 CONTINUE 
SIOP 

C ********JJ** ************************************************* ******* 

SUBROUTINE EtiE M (IELKM) 

COM MO M/Al /CORD (1C 1,2) , W>D 

CO tf MO N/A2/K0NECTC 100,2) f^ELEM 
COMMO.v/A3/ELGR(2 f 2),CLF(2) 

COMMO\/A4/DXX,R£TA,CST,XLAMD,0,UX,POP f DT 
CD w KON/AS/CONCClOil s „ „ 

COMMON/ A6/GLGK (1^1 , i Oi) . GLF( lyi ) 

COMMON / A’ 1 / NSPCD#IS p CO(2),SPCD(2) 

DIMENSION^ N0DEC2) , A ( ? , 2 > , ELGP ( 2 , 2 ) 

1 NODE (I ) =^0HECT(IELEM,1) 

H = 0 * 0 

DO 2 isl 7 

7 H = H?f COR^CHOntCn ,I)-CnRD(N0DE(2) 

H = 5DRr ( rO 
DO .3 1=1,2 
ELRCX)=O.D 
DO 3 J = i , 2 „ 

EOGKf X , J)=0.0 
ELGPCi, J)=G.O 

C 5o"ca2cULATE E0E*I£NTAI. ii MAJPXC|5 

ELGKC 1 , 1 )=DXX/(1 
ELGKCl ,2)=-DXX/C 
ELGK(2,l)=~0XX/( 

EL0K(2,2)=DXX/(1 Atn fi 

ELGP ( 1 , l)=CRETA*H)/( 3.0*DI) 

ELGP ( 1 , 2) = CRE:TA*H)/( 6.0*01) 

ELGP(2rl)=(RETA*H)/(6.0*0T) 

ELGP (2 r 2) = (RETA*H)/(3. 0*»x j 
ELF( 1 )=CO*CST*H)/(?.0*PQR) 

ELF(2)=(0*CST*H)/ (2.9*PQR) 





DO 1 lsi,NMQT> 

GLF(I)sGLFCI)-SPCO( J)*GLGK(-I,lSpCDf J) ) 

DO 3 llsl-NS^CO 

I=TSPCO(tI) 

GbF(I)=SPC n CTi) 

00 2 J=1,nN0D 

GLiGKCI, J)=0.0 

GLGKC J,I)=0.0 

GbGK(I,I)=1.0 

RETURN 

END 



uuuuuuuu 


5S55o55SS5S3RS?5«8ft;'55T55oS l,eM,? “ TS 

cB«3S^3^cg^(^ ST ' xlJA ’ ,t) ' a ' ux ' P3R ' ni 

CD«MON/Al/CQRi)({01 f 2)7N5ofl 12 
CDMH0N/A3/F°GK ( I!! ) 'Ib^C 3 ) EM 

cd«ho:;/a6/gl,gk t ioi , 101 ) , glfc 101 1 
CDM«aN/A7/NSPCO,ISPCO(?5,SPCO(2 

DATA 3 D x X,KFrA > c£x,XUHO ( U, y x, P GK,0 T/ 0. 3 ,l.C,0.0,0.0,r.O,n.n3u, 

DATA n£p?G,ISFCO/2. 1.101/ 

OATA NMon.t’iaKvm.So/ 

DATA SPCO/l(l.v,0.0/ 

DO 1 1=1. KNUD 

^^ c ?ii*^r u ' j) ' jsi ' 23 '^ rcT) 

READ (21 ,*) , (KQNECTCI, J) ,Jrl ,3) 

T = Q,"0 

DO 3 1=1,19 
T=T+DI 
CALL AMRb 
CALL SPCO.vC 

CALL, GAUSjn(CLGK f GLFVJ»*3P) 

TXPE *,NNOD 
DO 49 Jsl.NWnp 
C3NCC J)=f,Lr(J) 

WRiret47 ,*) ,t, J,canc( J) 

CONTINUE 

vSrnp 

END 

C* ******** **************************************************** ******* 

SUBROUTINE ELF M (IRLem) 

CON M On/ A l Z^ORlK m,2 l.NNUO 


49 

3 


common/ AJ/ r URi)U 01,2 },NNDD 
COMNDN/A2/KONECTC50,3) ,NELEM 
CDMHOw/A3/ELGKC3#3),ELF(3) 
C3MMnN/A4/nxy,RETArCSI,XLAMD,3#ur # PDP, 
C3MMDN/AS/CUNCC10J 5 



DT 


D r J 1 1 = 1. J 

NODECDaKCiMECTdELKM, I) 

H = 0.0 
DU 2 1=1.2 

H=N+( CORD ( NODE f 3 ) ,T)-C0RD(ND0E(3) ,1) )**2 
HsSQRrCrt) 

DO 3 1=1,3 
ELFCI)=0,0 


DO 3 J=l f i 
EL.GKCT, JJaO.O 


ELGPCI. J)=0.r 

CONTINUE ^ 

EbGKa?n = C7*Svx)/(3!o*H)-5x/(2.0*P3a) + CXLAi-iD*»ETA*n*2.0)/lb.E 
EbGK(i^2;=-( R,0*DX* )/(3.0*Hj+UX/(i * 5*P0R) + C XLAM£i*hEl A*H) /I % *> 
FLCKf 1 3^=[)/X/f 3.fr*/o-u£/(6.0*HDP)-lV|jAMn*PeTA*H) / 3( .0 , t 

c, r<( 2 # n--fR / *DXX) /<3#0*rf)-UX / (1 ♦ 5 *° 3 R> + (XLA'*» 0 *Pi:T**H)/lh.l 

PuGKf2 r 2?5flf. , y*DXX)/cS^O?Hji-(«.oixi l AviO*PETA*rt>/l5.v^ _ 

Pf ~K ( 2 * n — r p't n y y ) / ( 3 • 0 * ri ^ * U X / ( 3 • 5 * P 3 R ) + ( X L AA D * R t, T £ * f i ) / 1 S • L 
r b r<r V i ?-nyx/( 3 D*fn + UX/(f>!oiPoR)-CYuAHDipETA*H)/3t.o e 

‘5^555 ! # CtDxX)/C3.oiH)-UX/Cl,5^PnR) + CXLA M r^kFlA+U)/l5 i .f‘ 

lLGK(3',3) = fk'*DX$^f3.6*H)it'X/(7.0*PjR) + (2.0«X l .A«D*RETA» n l/l5. 

KbG?aai=CR | :i’A+H)/c7.5*Dn 
EDSPCI 21 = CEF.TA*H)/(lS.0*Dr) 

ELGPC 1 , 3js-(P£TA*H)/(30.0*DI) 

F.LGPC21 1) = CKETA*H)/(15.0*DT) 

ELCPC2, 2)=(3.u*RETA*H)/(15.p*pT) 

ELGP C 2 , 31=( RETAIN) /(1§*0*DT1 

ELG?C3,1)=*CPE t h*H)/(30.0*DX) 

ELGP(3,2)=(RETA*H)/( 15.0*pT) 

ELGPC3,3)=C2*RETA*H)/(15«0*DT) 

ELF (1 )a(0»CST*H)/C6.0*P0R) 

Et>C2)*U*CST*H)/(6.0*POR) 

ELF(3)=(D*CST*N)/(6,C*POR) 

XH=0.B75 
DD 5 1=1,3 

A? T 5 J5 = (f*LCPCl, J)-((1-TH)*ELGK(I, J) V 

ELG^CI, J) = (Et.GP(I, j) + rH*ELGKCI# j5) 

CONTINUE 
DO 6 1=1,3 



CDNriNul LPCIUA(T ' J)?cnf ' C(K0NECTflF:i,FM# J)) 

RETURN 
END 

* * * * * , T T T . r T T ^ ^ 
subroutine asmbl 
CD* MOiWAi/COR 0(101,21,^00 
CD*M0\I/A2/K0N£CT(50,3 5 .NELEM 
CDM*0 .\i/A3/ELGKC3 

C3MMON/A4/OAX. RETA,C5T,XE&Mr) f o # ity p 1R nr 
CDMM!DM/A5/CUn£(10i5 U, j,ux,P3R,nr 

CD V H0M/A6/GLGK(101 . 1 0 1 ) # GLFi 1 01 1 

C3MH0^/A7/NSPCD # JSPcnC2j»SPC0C7) 
do i i = i , I'iNun w 1 ' 

GLFCIJSO.O 

DO 1 jsi wwon 

GLGK(I. J)=9.0 

CONTINUE 

DD 2 l=X f NFLErt 

IEliEMsl 

CALL EL£M(TELE M ) 

DD 2 J=1 ,3 

g^CKDNECT(IELEM,J)) = GLF(KONECX c TELE^,J ) ) + ELF(J> 
CONTINUE 

return 

END 


5 

13 


9 

15 


DJ 1 1=1 | N 
BXGsA&S ( A C 1 » T ) ) ; JNil>BTG=l 
N 1 1 = M • I 1 ; T p ( T 1 • ijT . 2 ) GD Tj 3 
DD 2 J = 2 , N T 1 

IF(ARS(ACJ , X) ) •LE.HlG) GO TO 2 
BIG=AtfR(A( J,T) ) ; I«f>bIG=J 
CONTINUE 

IFClNDBin.FG.l) GO 10 3. 

T = C (I NDBTG) ;C(INDB 1 G)=CC 1 ) ;CC 1 I=T 
DD 3 K = 1 . N 

T = A(INDBIG#K) i A(INI>BTG,K)=A( 1 ,K)?aU,K)=T 
C J K T T N U E 
I P3 =1 + 1 

IFCIP1.3T.N) GO TO 13 
DD b KL=Tpl,N 
A ( J f KL)=AU ,KL)/A(1 ,1) 

C ( 1 )=CC1)/AC1,T) 

A(1,I)=1, 

DD 9 liMstf.N 

CCLM)=CCLA)-CC1 )*ACL",I) 

IFdPhST.ND GO TO 9 
IFCA(bn.l5,EO.0.U) GO TD 4 
DD 4 L=IPl#r* # , , 

ACLN,L)=A(LP ! ,L)-A(1 ,rO»A(L H ,I) 

CONTINUE 

ACLR,n=A(LK,n-ACi ,I)*A(L#,I) 

TC = C(1) 

DD 7 B\’ = 1,U 
T=A tl, w .O 

IF(MN.^I.I) C(MN-l)=C(MM 
DD b NOsJLu 
ACM3-1, ^U)=A(Un # MD) 

A(N BN)=T 

ctu)=rc 

CONTINUE 
RETURN 
END 

**♦ * V T f T T ^ T T T t ▼ 'r- 

subroutine: spconc 

C0Brt0N/Al/C0Rp(10l,23.«Ngp.. M 

C0KK0N/ft2/K0 , JEXT(5O, 3) , htUtH 
C3MM0,'./A3/EIiRKC3,3),Kl;F(3) • „ p R 

C0FH0N/A4/DXX.RE.TA, CnT,Xl'A f D»3rU*» p 3R;PX 

CDMHON/A5/CU'''CClOi) 

CORrtDN’/AfiXGLinKf 101 , 101 ) -GLf lIOl) 
caRHn>j/A7/«SPc - o,is[ , cn(2) ,s p co( 2 ) 

03 1 J=i,NSPCU 

Rl!r(I) = 3tF(I)-SPC0(J)*GbGKU,ISpCDCJ)) 

DO 3 Ilsl.NSPCO 
IsTSPCOCIi) 

GUF(I)=SPC0CT1) 

00 2 J*l,N«Op 
GLGKf 1, J)=0.0 
G3GKC J, I) = C ,0 



3 


GLGKCI 

RETURN 

end 


I)=l.o 



1 

2 

11 


49 

3 

12 


C 

C 

c 


NbOND= NUMBER OF' BOUNDARY nodes * * ' * ******* *** * ♦ 

^5S58;5W5}S4S2W!S.{5S4SBl! r »«*ge* m*d ^dary M3nE5! 

MATN PROGRAM FDk TWO DIMENSIONAL**!^** *»4^^»^»4f**»*4****»l*?** 
C0MMn,S/A!(/CDRD(147,2),NNU0 w 5 

COMMON/A?/KONECTC240,S) r N ELEM 

C0MM0N/A3/FbGK(j,3) , RLF( 3 ) ,LOCNnDf?').r>(A ■j ^ r / ^ 
CD«^DN/A5/CQNtU4^) VX ' DyV,fteT4 ' CKT ' X,,AMD ' Cl ^ X ’ 0 l'< r, ' 1H 'i' T ^^ P i 

C jMMOn / A 6/N fiflfiO. I Bnr;D(?0t 3) 

COMMON/ A7/G0Gn( I 47,147 ),& 0FC1 47) 

CO«HD^/ftP/NSPCO A lsfcO(lS),Spc5li4l 

r t) YX # DY Y . RET A # CST . XL A Mr . 0 , i?X , UY , PO D nT AN PH/ 

1 . 00 506, o.o.o.y,. 000009,1 .0,6. 0,0.0, 6.0 Ob. ol?"-td 1 c 1/ 
DATA «Nan,HE(.f«/147,240/ * * 1 ' " J ' 1 • »• 1 7 

DATAt^SPCU.I-SPGO/IA,! ,22, 4i, (>4,85,1 06,1 27, 21, 42,63, 

DATA NB5rJD/2f>/*' 5 *' 5 " 5 " S *' 5 ’' tj *'' J ' 0,U ' U * U ' e * 0 ' ( '*°'"*'"' J * l “ / 

03 1 1 = 1,I1A'0D 

READC22j*)(CnRD(I,J),J=i,2),C3i C fn 

si'iPid’isfe” 11 ''"'' 1 ' 1 ’” 

Tsfj.O 

DO 3 1=1, SO 
X = T + P X 
CALL A G M R L 
Call »Pcnrir 

C Alb GAUfi JD(GLGK.GLF,NNOD) 
type * , n k o n 

DO 49 Jsl.NUni) 

CDNCC J)=CLF(. J) 

WR7XFC43,*) O r COA'CfJ) 

COLXTMlIt: 

continue 

STOP 
KM D 

* #* * i 

SUBROUTINE bPCUhC 
CJM MON/hJ /CORDf I47 r 2) f «mud 
common /a?/ko,uct< 24 0,5) ,nelem 

C0 M MDN/A3/FbGK(3i’3) / LLF(.3),LnCNpDt2),0(3,3),F*t3) 

co M MON/A4/nxx,DXY,Dyy ,oyy ,reta,csi , xlamd,d, ux f uy , pd?.,dt, alph 
COMKON/AS/CONC f 147 J 
tOMMON/AS/NbOf’O. } BOLT) 1 2 0/3) 

COMMON / H 7/GLGK( 147, J47),&LF(147) 

COMMO t ‘i/AR/NSPCO,ISPCOCl^) ,SPCDtl4) 

Do 1 j=i,nspco 
DO 1 1=1 , LNOO 

GLFCI) = GLF(1 )-SPCO(0)*GLGK(If lSpraf J) J 
00 3 11=1 ,NSC>C0 
jsXbPcncTi) 

GLFtn = 3RC0Ui) 

03 2 J = 1 , »■<’ N 0 0 
GLGKfl, JJsC.C* 

GLGKC J,I)=r .0 
0LGNC1,I)=1 
RETURN 
END 

. . . 

SUbRi/UTIMf 5 'TkFaCC Tf’ODl , J ! 4*)DJ ) 

CDMH0%/Al/rijPD(147,2) f fW, 

C0 M MnA/A?/KONFCTf 240r J) # PC.LEJJ,, w .. p,., 

CO V MO i/m3/FLGK( 3, 3 ) ,f:LF( 3) ,LOCN3OC?)/0Ci,3) , 

C J M H0N/^4/nxy , L>x{ , of X , Dy v , RET A , COX , XLA*D# J, uX, L l ,POp. , pT , .- 
copmon/a^/cj m c:ci47) 

C0«M0.,/A6/L’B0iJ0 f lHO?iO(2Cf i) 

s;s?aKft»* }ws«»b«J 3 ') 


********** 


L p n 


x=ibonpci£lf* i ,2) 

: e< ‘ I® j i 

g 2'? =cnxxir cP/c§ir.Yx* /2.uH(u^r)-cuy*B))/c?.o^o^ 

p i5 = lAuP f «/2N-)*ccnxT*Cj-(Dyy«sn 

DO 21 Tl=i,J 

E(TI)=0.0 

DO 21 JJ=1 ,3 

DCII,J.U=0.l» 



21 CONTTNlIt; 

03 7 b=l , 2 
Tib-LDCM jficL) 

ECLb^KCIiD+PcL) 

03 7 M=l,2 
MMsbnCNDDCM) 

7 ^^i^r ctt " MM)+S(L ' , ' ) 

RETURN 
Cain 

SUBROUTING GAUSjOfA.C.N) 

OX«eNSIDW AU47 , 147 } f fc 147) 

Du 1 1=1 . N 

BISSABSCACI.IJISINDBICSI 

20 73 3 

J* J - 4 # N 1 1 

IFfABSf ACJ # X)),LE.8lG) GO TD 2 

2 S 88 Stt*«'»”»«ww 

iFriNDBic.EQ.i) go rn 3 
T*ClXHnBIG);C(INDBIC)=CCl)fC(l) S T 

P vJ 3 K~X ,h 

3 cjNrnll^ G,K) ' AtlNI)B T G ' K)sM1,K)?ACl,K ^ T 
TPIalfl - 

IFCIPl.GT.N) GO TO 13 
03 5 KIjsTP1,N 

5 A(1 ,KL)=A(l r KL)/A(l 7) 

13 CC 1 } = C(n/ACl f l) 

Ail # T ) = i * 

03 9 b w =2,N 

ccr»^)=r(LM)-rci)»ACLM r i) 

IFCJOl.^.N) GO TO 9 
IF(bCLM, D.EO.O.O) GO TO 4 
DO 4 L = X P 1 * K 

ACCjM.uIsACIiM, D-AC1 ,I.)*A(L'M) 

4 CJ^TI^^E 

P A t , 1 ) “A Ui'W 1 ) -A ( 1 , J ) * A C L i4 ,X ) 

is rcsecn 

DO 7 MU=t,N 
r=A ti f nt) 

IF(MN.Gt.X) C(MN-1)=C(MM) 

03 b Nn = ? f fJ 

6 ACMO-l f MM)=ACNO # MN) 

7 A(U, w iOsT. 

cm=rr 

1 CONTINUE 

RETURN 
EnIO 

4 4 * 1 

SUBROUTINE fciM£«(lEliEPj 
C COMMON/*! /C0PD(l»N0D-2) , H N QD 

C^HMnw/A\/CDRD(J47,J) ,MNOD 
C C3M nD U / A2/K0NECT C RELEfi f 3 ) # fJ ET'E M * . 

CU^ , rtO\/A2/KO' T ECT(^4G JL 3),u£LEM 

C3WKO\/A3/EJ.GKC3#3) , FbF( 3) , L0C*3D ( ?) , D f 3 ,3 5 , F( 3) 

C3*>lON/A4/DXX # DXY,DYX,DyY,KETA,CST,XLA’ , D,Q,UX,ny ,pnk,,.T,ALPH 

C CD^MDN/AS/CDN’CCNNOO) 

CjM*\0v/AS/C0NCC147) 

C COLSON 'Afi/NbOND, IbOWDC NBOfcD, 3) 

CO^ v tON/HG/NbOND £ lBONOC?Or3) 

01 MEN SIjn N()DE(3) ,A(3 r 3) ,EbGP(i,3) 

1 NDOtjrnsftON'ECTdEbEMf X) 

Bl=C0ROCW00Ef2j f 2)-cr|R0(ND , >EC3K2) 

B2 = COP.O(OOOE(3) ,?)-COROCNuOE(i) ,?) 

BJaCOROCNOOEf 1 ) , 2 ) -CORO ( N ODE ( 2 ) , ? ) 

Cl=CDf:0(NU0EC3) , 1 ) -CORO (NODE (2) ,1) 
cJ-COKOCNOOEC 1) I 1 )-COrO(NOO|C 33 r n 
C^C0KOt f J0OEC2) f l)"COKDCN3DECl) ,1 j 

EA»( CC1 *O3)-(C3*Bl))/2.0 
ea=absc ea) 

DO 3 1 = 1,3 
EGECIlsaiO 
03 3 jsl ,i A 
EbGK ( X * J ; =0, 0 

EuGPU, J)s0.0 
3 CONTINUE 

C HAriacEs 



S 7 t o < X > r hy X * R { * ft? U 7 DX Y *B 1 * C? ) + C 0 Y **~1 *&2 ) + (. n * I 1 ¥ A ) J / 
lC^oi£A);((UX^li+CW^Cn)/C£>.0*Pnbt(XLAVD*REXMCM/12.U 



<r>*y*C2**?3 )/ 


f f 1 1,1 ) = ( ( nXX*B3»ui i + f i)vYiH^l> l ? pns 5 + (XLiAM0*k' : "T! l ^FA 1 / I o 

pfrjVr^ + JJH x * bJ ) + » W3))/(fi Sl^^^^iJ^iyScsirin/ 

■ liiJiSjifscs 

I iiiiill ll 

ELGPCjJnsCHETAjEAJ/Cli’ulD?) 

£P£^ = < 0 * CST *EA)/C3;5*PORr 

£;!**£ £ ? > 55 ^ ) / { 3 m G*POR 3 

S4£i}^ = (^*CST* KA )/( 3 . 0* PCJR 3 
run^n ^O^NOAHy ELEMENT 

XgnMO=T3AR CH Cl2) L b b T 


r«n*i=iRow+i 

?2 2|>J I2=trdwi,wbomd 

I w“LlN2U 

s ° t:i 2qi 

I«n-t=x2 
gj rn ?u? 

CDwriMig 

cjMridne 

XFCirHECK.FO.iO Go TO 3i)l 

J2???.fI 30Hn < r R 0 W'2J ? JNOni S T B ONn<TBnj -o 
Iftxnani.M.rjNECTf m.rS.oT'PsWU??*'#* 


X^^JNjDi*E0.K0^ECrfIEtjEB#3)j uOr ViD ( 7 1 — ^ 

clit s6»F-grc?SAi£j:3Sgg?$ 2> ' iM - D i'J K 3Di _ 

G J ID 4 0» i 
CONTINUE 

Trial ). 9 * 

OJ b 1= 1,3 

DD S J=Ii3 

&?.rl } =y If J’ + cf.I'GPci , jo - u i . 0-rnjiEi.GK c t , j» j 

coririNi'E 5 au,PtI,J)+ * H tElj3KCJ,J ^ )) 

2? t. 1 = 1,3 

no 6 J = i,3 

E^^^pCn + RiiFCJJ+AfJL, Jj^CDMCCKONFCTaELEM^i)} 

rltur* r ' 

EtiO 

SUnRnjTIJJE *^*\p^* ***************************** ************ ******** 

CO^ON/Ai/ruRijC 147,2) ,WNDD 

rO N IHO\/A?/KO«fcCT(24o f S) # NEIiEM 

£u*J ■/A3 / rn r ,l';C3,3) / F:t,F"C3) ,LOCNoP(7) ,p( 3,3) ,F(3) 

C323S^J^CO^f»4f) fX ' nn ' BEn ' CSI ' XUMD ' 0 ' UX ' nX ' P?K ' , ‘ T ' A, ' Prt 

CO'^n.^/rtT/GL^K (14 7,147) ,Gl,F(147) 
r?F?r4-, 1,, i Mon 

GL»f C T 3 = t».«> 
n u i jai-HNon 
GLGK ( 1 , J $ SO. G 
CSHfXNUB 
00 2 l = i,UELEl! 

ISI.E>*si 

CAOO ET.EK £ Tt.T.E;'* J 
00 2 J=1 , 5 

GOFCK JMEOXCIEI.tlM, J))=G0F c KOHRcr(ItU.EK,J))+El,F(J) 

GURKfKnrtfccT(XEr.Ef,»l),KnwBCX(IEUE»',K))sGbGK(KOAFCTCTEOE«l,.»), 

X Kn.’lECT C I lil'E" , K ) J + ELPKC J , K ) 
com n mis 
DO 71 1*1, h 

WKTX'E(2 0,*1 , C<;l-GKU,J),J=1,6) 

REFORM 

PWD 









Distance (metres) 




Distance (metres) 





Fig. 4.6 Comparison of analytical results and F.E.M. results for quadratic ca 




IMPERVIOUS 

FIG. 4-7(a) Concentration contours for two dimensional case at t = 10 mins. 
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